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ABSTRACT 

This paper describes the implementation of optimization techniques 
based on control theory for wing and wing-body design of supersonic 
configurations. The work represents an extension of our earlier re- 
search in which control theory is used to devise a design procedure 
that significantly reduces the computational cost by employing an 
adjoint equation. In previous studies it was shown that control 
theory could be used to devise transonic design methods for air- 
foils and wings in which the shape and the surrounding body-fitted 
mesh are both generated analytically, and the control is the map- 
ping function [5, 6, 8]. The method has also been implemented 
for both transonic potential flows and transonic flows governed by 
the Euler equations using an alternative formulation which employs 
numerically generated grids, so that it can treat more general config- 
urations [16, 9, 17]. Here results are presented for three-dimensional 
design cases subject to supersonic flows governed by the Euler equa- 
tion. 

INTRODUCTION 

Since the inception of CFD, researchers have sought not only accu- 
rate aerodynamic prediction methods for given configurations, but 
also design methods capable of creating new optimum configura- 
tions. Yet, while flow analysis can now be carried out over quite 
complex configurations using the Navier-Stokes equations with a 
high degree of confidence, direct CFD based design is still limited to 
simple two-dimensional and three-dimensional configurations, usu- 
ally without the inclusion of viscous effects. 

One approach to extend the existing CFD analysis capability to 
treat the design problem, but only at the price of heavy compu- 
tational expense, is to use numerical optimization methods. The 
essence of these methods is very simple: a numerical optimization 
procedure is coupled directly to an existing CFD analysis algorithm. 
The numerical optimization procedure attempts to extremize a cho- 
sen aerodynamic measure of merit which is evaluated by the given 
CFD code. The configuration is systematically modified through 
user specified design variables. Most of these optimization proce- 
dures require the gradient of the the objective function with respect 
to changes in the design variables. The simplest method of obtaining 


gradient information is by finite differences. In this technique, the 
gradient components are estimated by independently perturbing each 
design variable with a finite step, calculating the corresponding value 
of the objective function using CFD analysis, and forming the ratio 
of the differences. The gradient is used by the numerical optimiza- 
tion algorithm to calculate a search direction using steepest descent, 
conjugate gradient, or quasi-Newton techniques. After finding the 
minimum or maximum of the objective function along the search 
direction, the entire process is repeated until the gradient approaches 
zero and further improvement is impossible. 

The use of numerical optimization for transonic aerodynamic 
shape design was pioneered by Hicks, Murman and Vanderplaats 
[4]. They applied the method to two-dimensional profile design sub- 
ject to the potential flow equation. The method was quickly extended 
to wing design by Hicks and Henne [3], More recently, in the work 
of Reuther, Cliff, Hicks and Van Dam, the method has proven to be 
successful for the design of supersonic wing-body transport config- 
urations [15]. In all of these cases finite difference methods were 
used to obtain the required gradient information. 

FORMULATION OF THE DESIGN PROBLEM AS A CON- 
TROL PROBLEM 

In an attempt to reduce the computational expense of aerodynamic 
optimization recent works have focused on using control theory to 
provide inexpensive gradient information. The technique can be 
outlined as follows. 

For flow about an airfoil or wing, the aerodynamic properties 
which define the cost function are functions of the flow-field vari- 
ables (w) and the physical location of the boundary, which may be 
represented by the function T, say. Then 

I = l(w,D 

and a change in results in a change 

in the cost function. Each term in (1), except for 6w, can be easily 
obtained. Finite difference methods require a number of additional 
flow calculations equal to the number of design variables. Using 
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control theory, the governing equations of the flowfield are intro- 
duced as a constraint in such a way that the final evaluation of the 
gradient does not require multiple flow solutions. In order to achieve 
this, 6w must be eliminated fiom (1). The governing equation R 
and its first variation express the dependence of w and T within the 
flowfield domain D: 


R(w,:F) = 0, 



dT 


8T = 0. 


( 2 ) 


Next, introducing a Lagrange multiplier we have 



Choosing to satisfy the adjoint equation 



the first term is eliminated, and we find that the desired gradient is 
given by 



Since (4) is independent of Sw, the gradient of 1 with respect to an 
arbitrary number of design variables can be determined without the 
need for additional flow-field evaluations. The main cost is in solving 
the adjoint equation (3). In general, the adjoint problem is about as 
complex as a flow solution. If the number of design variables is 
large, the cost differential between one adjoint solution and the large 
number of flowfield evaluations required to determine the gradient 
by finite differencing becomes compelling. Once equation (4) is 
obtained, Q can be fed into any numerical optimization algorithm to 
obtain an improved design. 


ISSUES OF IMPORTANCE FOR DESIGN PROBLEMS 

The development of aerodynamic design procedures that employ an 
adjoint equation formulation is currently being investigated by many 
researchers. These methods promise to allow computational fluid 
dynamics methods to become true aerodynamic design methods. 
References [1, 2, 10, 13, 12, 11, 14, 19] represent a partial list of 
recent works that reflect this developing field. However, as is the 
case in any new research field, many questions remain. Probably the 
most salient issues of concern are the following: 

1. Discrete vs. continuous sensitivities 

2. Choice of optimization procedures 

3. Treatment of geometric and aerodynamic constraints 

4. The level of coupling between design and analysis 

5. The parameterization of the design space 

These topics still require intensive investigation in the up coming 
years. With regards to the first item, it is historically interesting that 
Jameson in 1988 [5], first developed the equations necessary for a 
continuous sensitivity approach to treat the design of airfoils and 


wings subject to transonic flows. By continuous sensitivities it is 
implied that the steps represented by equations (1-4) are applied to 
the governing differential equations. The resulting adjoint differen- 
tial equations with the appropriate boundary conditions may then be 
discretized and solved in a manner similar to that used for the flow 
solver. One may alternatively derive a set of discrete adjoint equa- 
tions directly from the discrete approximation to the flow equations 
by following the procedure outlined in equations (1-4). The resulting 
discrete adjoint equations are one of the possible discretizations of 
the continuous adjoint equations. This alternative is mentioned in 
[6], but was not adopted in that work because of the complexity of 
the resulting discrete adjoint system. The approach has found favor 
in the work of Taylor et al. [13, 10] and Baysal et al. [1, 2]. 

It seems that both alternatives have some advantages. The con- 
tinuous approach gives the researcher some hope for an intuitive 
understanding of the adjoint system and its related boundary condi- 
tions. The discrete approach, in theory, maintains perfect algebredc 
consistency at the discrete level. If properly implemented it will 
give gradients which closely match those obtained through finite dif- 
ferences. The continuous formulation produces slightly inaccurate 
gradients due to differences in the discretization. However, these 
inaccuracies must vanish as the mesh width is reduced. 

The discrete adjoint equations are a linear system, whether derived 
directly by the discrete approach or by discretizaton of the continuous 
adjoint equation. The size and complexity of the system, however, 
makes the use of direct solution methods unrealistic for all but the 
smallest problems. The application of the continuous sensitivity 
analysis fosters the easy recycling of the flow solution algorithm fot 
the solution of the adjoint equations, since the steps applied to the 
original governing differential equations can be duplicated for the 
adjoint differential equations. When the discrete approach is used, 
the adjoint equations have a complexity which makes it hard to find 
decompositions to facilitate their solution unless the structure from 
the continuous adjoint is used as a guide. The discrete method is 
subject, moreover, to the difficulty that the discrete flow equations 
often contain nonlinear flux limiting functions which are not differ- 
entiable. It also limits the flexibility to use adaptive discretization 
techniques with order and mesh refinement, such as the /i —p method, 
because the adjoint discretization is fixed by the flow discretization. 

It turns out that the determination of items (2-4) in the above fist 
strongly hinge on the choice for (5). In Jameson’s first works in the 
area [5, 6, 8], every surface mesh point was used as a design variable. 
In three-dimensional wing design cases this led to as many as 4224 
design variables [8]. The use of the adjoint method eliminated 
the unacceptable costs that such a large number of design variables 
would incur for traditional finite difference methods. If the approach 
were extended to treat complete aircraft configurations, at least tens 
of thousands of design variables would be necessary. Such a large 
number of design variables precludes the use of descent algorithms 
such as Newton or quasi-Newton approaches simply because of the 
high cost of matrix operations for such methods. The limitation of 
using a simple descent procedure, such as steepest descent, has the 
consequence that significant errors can be tolerated initially in the 
gradient evaluation. Such methods therefore favor tighter coupling 
of the flow solver, the adjoint solver, and the overall design problem 
to accelerate convergence. Ta’asan et al. [11], have taken advantage 
of this by formulating the design problem as a one shot procedure 
where all three systems are advanced simultaneously. Choosing such 
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a design space suffers from admitting poorly conditioned design 
problems. This is best exemplified by the case where only one point 
on the surface of an airfoil is moved, resulting in a highly nonlinear 
design response. In his original work Jameson [ 6 , 7 ] addressed this 
poor conditioning of the design problem by smoothing the control 
(surface shape) and thus removing the problematic high frequency 
content fix>m the advancing solution shape. 

Hicks and others [ 4 , 3 , 15 ] have in the past parameterized the 
design space using sets of smooth functions that perturb the initial 
geometry. By using such a parameterization it is possible to work 
with considerably fewer design variables than the choice of every 
mesh point And since Hick’s original works exclusively used finite 
difference gradients, the inherent limit to a small number of design 
variables allowed for the use of more efficient search strategies. 

When distributed over the entire chord on both upper and lower 
surfaces, such analytic perturbation functions admit a large possible 
design space. They can be chosen such that symmetry, thickness, 
or volume can be explicitly constrained, thus avoiding the use of 
expensive constrained optimization algorithms to address geometric 
constraints. However, such a choice of design variables does not 
guarantee that a solution, for example, of the inverse problem for 
a realizable target pressure distribution will necessarily be attained. 
Finally, they have one last advantage over using the mesh points; 
there is no need to smooth the resulting solutions as the design 
proceeds when using Hicks-Henne functions, since by construction 
higher frequencies are not admitted and thus the design spaces are 
naturally well posed. 


Introduce contravariant velocity components as 



The Euler equations can now be written as 


with 


dW dFj 
dt d^i 


0 in 


( 9 ) 
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Assume now that designs are limited to wing and wing-body con- 
figurations using a C-H topology mesh. If a body is present it is 
assumed to exist as part of the face containing the symmetry plane. 
Thus the new computational coordinate system conforms to the wing 
in such a way that the wing surface Bw is represented by ^2 = 0 
and the body and symmetry plane Bb conform such that fs = 0 . 
Then the flow is determined as the steady state solution of equation 
( 9 ) subject to the flow tangency conditions 


THREE-DIMENSIONAL DESIGN USING THE EULER 
EQUATIONS 


The application of control theory to aerodynamic design problems 
is illustrated by treating the case of three-dimensional wing design, 
using the inviscid Euler equations as the mathematical model for 
compressible flow. In this case it proves convenient to denote the 
Cartesian coordinates and velocity components by xi , X2, 13 and tti , 
U3* The three-dimensional Euler equations may be written as 


where 


— , 

dt dxi 


0 in 
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' pUi 

pu\ 


puiU\ -\-p8ii 

PU2 

> , fi = < 

pUiU2 -h p8i2 

pm 


puim +p^i3 

. P^ . 




and Sij is the Kronecker delta function. Also, 

p = (7 («?)}, 


( 5 ) 


( 6 ) 


( 7 ) 


and 

pH = pE+p (8) 

where 7 is the ratio of the specific heats. Consider a transformation 
to coordinates fi, (2. 6 where 





J = det(if), A'-’ 


dxj 


1/2=0 onBw, C/ 3=0 onBs. (11) 

At the far field boundary Bpy freestream conditions are specified 
for incoming waves, while outgoing waves are determined by the 
solution. 

In our previous works, the illustrative problem most often used 
specified the cost function as the difference between the current and 
some desired pressure distribution. However for supersonic flows 
where pressure distributions for optimum aerodynamic performance 
are less clearly understood then for transonic flows, it is advantageous 
to consider the use of drag as the cost function. 

I = Cdw + 

= {Caw + Cab) cos a + {Cnw + Cnb ) « 

= ^ I I Cp ( 5 x cos a + sin O') 

r©i 

+ — f f Cp cos a Sy sin a) d^\d^2y 

^ref J Jbb 

where Sx and Sy define projected surface areas, and is the 
reference area. The design problem is now treated as a control 
problem where the control function is the wing shape, which is to be 
chosen to minimize I subject to the constraints defined by the flow 
equations ( 5 - 10 ). A variation in the shape will cause a variation 6 p 
in the pressure and consequently a variation in the cost function 

61 = 6Cy^ cos a -I- 6Cs sin a 
-f {-Ca sin a -I- Cjq cos a} 6a 

fdCA . 1 , 

-I- < cos Q + sin a ^ 6a. 

L 5a da J 
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where 6Ca and 8Cf^ the variation due to changes in the design 
parameters with a fixed. To treat the interesting problem of prac- 
tical supersonic design wave drag must be minimized at a fixed lift 
coefficient Thus an additional constraint is given by 


which gives 


+ 

+ 


8Ci, = 0 = + bCtf. 


SCj^ cos a — SC A sin o 
\-Cff sin a - Ca cos a) Set 


(££iLcosa-^sma|5a=0 
I 9a 9a J 

Combining these two expressions to eliminate Sa gives 

61 = SC A cos a -I- 6C]^ sin a 

-l-li {SCn cos a - SC A sin a} , 

where Q is given by 

(-CL + ^cosg+^sina) 
(Cu + ^cosa-^sina) ‘ 


( 12 ) 


Since p depends on to through the equation of state (7-8), the varia- 
tion Sp can be determined from the variation If a fixed computa- 

tional domain is used the variations in the shape result in variations 
in the mapping derivatives. Define the Jacobian matrices 


A = |J. Ci = JK-/Ai. (13) 

Then the equation for in the steady state becomes 

4 - = 0 . 

where 

SFi=CiSw + s(^J^"^ fj. 

Now, multiplying by a vector co-state variable assuming the result 
is differentiable, and integrating by parts over the domain 

where n, are components of a unit vector normal to the boundary. 
The variation in the cost function can also be expressed in terms of 
Sp after (12 and 14) are summed to give, 


+ — If Cp{ (55a: cos a -I- SSy sin a) 

^ref JJbs 

-|-D (55y cos a — SSx sin a^ | d^\ 

- J + J 

On the wing surface Bwy ni = ns = 0 and on the body and 
symmetry plane 5 b, ni = n 2 = 0. It follows from equation (11) 
that 


SF2 = J < 


55s = J < 


0 


0 

l^^p 


>+p< 




«(^e) 


‘(■'ij) 

‘(■'ft) 


> on Bw 


> on Bb‘ 


(16) 


Suppose now that ^ is the steady state solution of the adjoint 
equation 

in D. (17) 


^-cT^=o 

o/- 


dt ’ d^i 

For supersonic flows the choice of boundary conditions at the outer 
domain can be developed from physical intuition as well as mathe- 
matics. For a given geometry, say a wing, a change in the surface at 
any particular point, V, will incur changes in the flowfield and hence 
the performance in the area defined by the Mach cone originating at 
“P. Similarly, it is possible to determine the region over which sur- 
face changes effect the flow condition at a given point This region 
would also form a cone that would point roughly in the opposite 
direction, depending on local conditions, as the Mach cone. It is 
the solution of this reverse problem that the adjoint represents. The 
contribution to say drag, at a given point is influenced by changes to 
the surface at all points within the reverse cone. The correct farfield 
boundary conditions for the adjoint equation that are consistent with 
this reversed character are. 


61 = 


COS a — 5x sin a)}d<,d6 

+ i // fip { (5* cos o + 5„ sin a) 

(s, cos a — Sx sin o) | d (3 

-I- — / / C7p I (55x COS a -I- 55y sin a) 

^ref J Jbw 

-l-D {SSy cos a - SSx sin a) | di\ diz 


^i _5 = 0 at the exit 

i!f\ -5 extrapolated from the interior at the entrance 


Then if the coordinate transformation is such that S{^JK * ) is neg- 
ligible in the far field, the last integral in (15) reduces to 



il)^SF2 d(id^3 



t/f'^SFz d^id^2- 


Thus by letting rp satisfy the boundary conditions. 


( 18 ) 


J 


f , ^^2 , di2 , 


^>1 

9X3/ 


Q on 5 IV and^B, 


(19) 
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where 


1 

+Q (5y COS a - 5x sin a) } 

We find after integrating by parts again that 

SI — — f j C^p / (SSx cos CK ^ sin of^ 

JjByy 

cos a — sin I d^3 

+ J J Cp { (SSx cos a + SSy sin a) 

+Q (55y cos a - 5S* sin a) } di^ ci^2 
+ j ^'^■^,{SGijSj)diu ( 20 ) 

where Gi> = 

In our past works it was shown that the remaining variational 
terms in (20) that are related to mesh and surface deformation, can 
be obtained through either an analytic mapping procedure [6, 8] or an 
analytic mesh perturbation procedure [16, 9, 17]. The entire design 
procedure is outlined as follows: 

1 . Solve the flowfield governing equations (5-10). 

2. Solve the adjoint equations (17) subject to the boundary condi- 
tion (19). 

3. Determine the mesh variational terms from analytic mesh per- 
turbation relationships. 

4. Constract SI for each design variable to form the gradient by 

( 20 ). 

5. Calculate the search direction based on a quasi-Newton algo- 
rithm and perform a line search. 

6. Return to (1). 

The basic method here builds on that used in reference [17] with 
the proper extensions to treat supersonic flow cases and more in- 
volved objective functions. Both the flow and adjoint solvers use an 
explicit multi-stage Runga-Kutta like algorithm accelerated by resid- 
ual smoothing and multigridding. Design variables were chosen as 
a set of Hicks-Henne functions that tend to eliminate the presense of 
high frequency response in the design space [18]. 

NUMERICAL TEST CASES AND CONCLUSIONS 

Two cases of design in the presense of supersonic flow are exam- 
ined. In the first case a wing alone design is carried out on a cranked 
delta plan form shown in Figure (1), with the body removed. The 
initial wing is defined as a flat, untwisted, and uncambered wing 
with a thickness distribution varying from 4.25 % at the root to 2.40 
% at the tip. The outboard wing panel is swept ahead of the Mach 
cone at the cruise Mach number of 2.4 and thus has a sharp leading 
edge. The inboard subsonic portion of the wing retains a blunt lead- 
ing edge to maintain low speed performance. The design variables 


{ (^Sx cos a -f- 5y sin a) 



1a: Wing-Body configuration 


Figure 1 : C-H Mesh Topology 193x41x49. 


consisted of 90 Hicks-Henne functions and 9 twist variables placed 
over the entire span and in such a way as to keep thickness constant 
The design problem involved minimizing total inviscid drag while 
keeping Cl constant at a value of 0.12. The flow solver was run 
in the constant lift mode by using an outer iteration to modify a. 
The necessary and terms used in the construction of the 
adjoint boundary condition are estimated by running two initial flow 
solutions at slightly different values of a. Figure (2) shows both the 
original and final airfoil sections and pressures at four representative 
stations. It is seen from these figures that the loading has shifted 
inboard while the shock on the leading edge inboard region ha been 
eliminated. There is also a trend to reduce the overall levels of min- 
imum pressure across the entire upper surface. The overall inviscid 
drag was reduced by 6.33% in 25 iterations, from Cd = 0.01121 to 
0.01050. 

In the second design example the complete wing-body configu- 
ration illustrated in Figure (1) was optimized. The configuration 
started with the same flat, uncambered, untwisted wing as that used 
in the wing alone design. In this case 121 design variables were 
used, consisting of 109 Hicks-Henne functions and 12 twist vari- 
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ables distributed over the wing. The fuselage was fixed during the 
design process. Lift was once agin fixed at Cl = .12 by iterating 
on alpha. Figure (3) shows both the original and final airfoil sec- 
tions and pressures at four stations. The overall inviscid drag was 
reduced by 6.07% in 8 iterations, from Cd = 0.01070 to 0.01005. 
The figures show that the same trend to reduce the minimum upper 
surface pressure is again seen for this wing-body design. The strong 
shock seen at the most inboard station of the initial design has been 
completely eliminated by the final design. Unlike the wing alone 
case the loading distribution is not drastically effected. Finally a 
small upper surface low pressure spike is retained in the final design 
at the leading edge and over the blunt leading edge. This spike has 
been seen in other design attempts of similar planforms indicating 
the possible advantage of a leading edge thrust component 

These results are comparable to those that could be obtained by 
previous methods using finite differences to calculate the gradients 

[15], and it appears that there is little room for further improvement 
at this flight condition with the given planform. 

In this research the development of an aerodynamic shape opti- 
mization method has been extended to treat the design of wing and 
wing-body configurations subject to supersonic flows. The method 
is roughly 30 times more efficient then comparable finite difference 
methods due to the use of an adjoint solver to obtain the gradi- 
ents. The clear limitation in these current results is the reliance on 
a single structured block for both the state and co-state fields. In 
the future our group intends not only to extend the method to treat 
both multiblock and unstructured meshes, but also to implement a 
Navier-Stokes version. 
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2c: span station z = 0.609 2d: span station z = 0.828 


Figure 2: SYN87 Wing Alone Design. 

M = 2.4, Cl = 0.120, InHial Cd = 0.01121, Final Cd = 0.01050 
90 HickS'Henne variables, 9 twist variables. 

— , o InHial Wing 
— , + Design after 25 iterations. 
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3a: span station z = 0.252 


5 1 



3b: span station z = 0.441 



3c: span station z = 0.653 


3d: span station z = 0.842 


Figure 3: SYN87 Wing-Body Design. 

M = 2.4, Cl = 0.120, Initial Co = 0.01070, Final Cd = 0.01005 
109 Hicks-Henne variables, 12 twist variables. 

— , o Initial Wing 
— , + Design after 8 iterations. 
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